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Abstract. The paper is devoted to the spectral analysis of effective preconditioners for linear 
systems obtained via a Finite Element approximation to diffusion-dominated convection-diffusion 
equations. We consider a model setting in which the structured finite element partition is made 
by equi-lateral triangles. Under such assumptions, if the problem is coercive, and the diffusive and 
convective coefficients are regular enough, then the proposed preconditioned matrix sequences exhibit 

■ a strong clustering at unity, the preconditioning matrix sequence and the original matrix sequence are 

■ spectrally equivalent, and the eigenvector matrices have a mild conditioning. The obtained results 
allow to show the optimality of the related preconditioned Krylov methods. The interest of such a 
study relies on the observation that automatic grid generators tend to construct equi-lateral triangles 
when the mesh is fine enough. Numerical tests, both on the model setting and in the non-structured 
case, show the effectiveness of the proposal and the correctness of the theoretical findings. 
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1. Introduction. The paper is concerned with the spectral and computational 
analysis of effective preconditioners for linear systems arising from Finite Element 
approximations to the elliptic convection-diffusion problem 



= 0, 



div ^-<z(x)Vit + /3(x)ttJ = /, x G 0, ^ ^ 

u\aa 



with SI domain of R 2 . We consider a model setting in which the structured finite 
element partition is made by equi-lateral triangles. The interest of such a partition 
relies on the observation that automatic grid generators tend to construct equi-lateral 
triangles when the mesh is fine enough. 

The analysis is performed having in mind two popular preconditioned Krylov 
methods. More precisely, we analyze the performances of the Preconditioned Con- 
jugate Gradient (PCG) method in the case of the diffusion problem and of the Pre- 
conditioned Generalized Minimal Residual (PGMRES) in the case of the convection- 
diffusion problem. 

We define the preconditioner as a combination of a basic (projected) Toeplitz 
matrix times diagonal structures. The diagonal part takes into account the variable 
coefficients in the operator of (1.1), and especially the diffusion coefficient a(x), while 
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the (projected) Tocplitz part derives from a special approximation of (1.1) when set- 
ting the diffusion coefficient to 1 and the convective velocity field to 0. Under such 
assumptions, if the problem is coercive, and the diffusive and convective coefficients 
are regular enough, then the proposed preconditioned matrix sequences have a strong 
clustering at unity, the preconditioning matrix sequence and the original matrix se- 
quence are spectrally equivalent, and the eigenvector matrices have a mild condition- 
ing. The obtained results allow to show the optimality of the related preconditioned 
Krylov methods. It is important to stress that interest of such a study relies on the 
observation that automatic grid generators tend to construct equi-lateral triangles 
when the mesh is fine enough. Numerical tests, both on the model setting and in the 
non-structured case, show the effectiveness of the proposal and the correctness of the 
theoretical findings. 

The outline of the paper is as follows. In Section 2 we report a brief description of 
the FE approximation of convection-diffusion equations and the preconditioncr defi- 
nition. Section 3 is devoted to the spectral analysis of the underlying preconditioned 
matrix sequences, in the case of structured uniform meshes. In Section 4, after a 
preliminary discussion on complexity issues, selected numerical tests illustrate the 
convergence properties stated in the former section and their extension under weak- 
ened assumption or in the case of unstructured meshes. A final Section 5 deals with 
perspectives and future works. 

2. Finite Element approximation and Preconditioning Strategy. Prob- 
lem (1.1) can be stated in variational form as follows: 



w here H%(Q,) is the space of square integrable functions, with square integrable weak 
derivatives vanishing on dfl. We assume that fl is a polygonal domain and we make 
the following hypotheses on the coefficients: 



The previous assumptions guarantee existence and uniqueness for problem (2.1) and 
hence the existence and uniqueness of the (weak) solution for problem (1.1). 

For the sake of simplicity, we restrict ourselves to linear finite element approxi- 
mation of problem (2.1). To this end, let Tk = {K} be a usual finite element partition 
of fl into triangles, with Kk = diam(if ) and h — maxjf Kk- Let V/, C Hq(Q) be the 
space of linear finite elements, i.e. 

Vh = {iph : f2 — > K s.t. (fh is continuous, <Ph\ K is linear, and <Ph\dn = 0}- 
The finite element approximation of problem (2.1) reads: 




(2.1) 




(fl), with a(x) > a > 0, 
(fl), with div(/3) > pointwise in fl 
(fl). 



(2.2) 



find Uh € Vh such that 

In [ a ^ u h ■ V<^ - P ■ Vip h u h ) = J n ffh for all ip h e V h . 



(2.3) 



For each internal node i of the mesh Th, let ipt € Vh be such that ^i(node i) = 1, 
and <^i(node j) — if i ^ j. Then, the collection of all </Vs is a base for Vh- We 
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will denote by n(h) the number of the internal nodes of Th, which corresponds to the 
dimension of Vh- Then, we write Uh as Uh = Y^j=i u j i Pj an d the variational equation 
(2.3) becomes an algebraic linear system: 

' fa, i = l,...,n(h). (2.4) 



y~] ( / • v ^ - $ ■ ^<Pi ) u j = / . 

~[ \Jn / Jn 



According to these notations and definitions, the algebraic equations in (2.4) can be 
rewritten in matrix form as the linear system 

A„(a,/3)x = b, A n {aJ) = 6„(a) + f„(/3) el" x ", n = n(h), (2.5) 

where 9„(a) and ^ n {P) represent the approximation of the diffusive term and ap- 
proximation of the convective term, respectively. More precisely, we have 

{@n(a))i,j = [ a VwVw, (2.6) 
Jn 

{9 n 0))ij =- Itf-Vtpi) <fij, (2.7) 
Jn 

where suitable quadrature formula are considered in the case of non constant coeffi- 
cient functions a and /3. 

As well known, the main drawback in the linear system resolution is due to the 
asymptotical ill-conditioning (i.e. very large for large dimensions), so that precon- 
ditioning is highly recommended. Hereafter, we refer to a preconditioning strategy 
previously analyzed in the case of FD/FE approximations of the diffusion problem 
[12, 15, 16, f 8, 19, 17] and recently applied to FD/FE approximations [3, 4, 10] of (f .1) 
with respect to the Preconditioned Hermitian and Skew-Hermitian Splitting (PHSS) 
method [2, 3]. More precisely, the preconditioning matrix sequence {P n (a)} is defined 
as 

P n (a) = Dl(a)A n (l,0)Dl(a) (2.8) 

where D n (a) = diag(A„(a, 0))diag _1 (A„(l, 0)), i.e., the suitable scaled main diagonal 
of A n (a, 0) and clearly A n (a, 0) equals Q n (a). 

The computational aspects of this preconditioning strategy with respect to Krylov 
methods will be discussed later in section 4.1. Here, preliminarily we want to stress as 
the preconditioner is tuned only with respect to the diffusion matrix 8 n (o): in other 
words, we are implicity assuming that the convection phenomenon is not dominant, 
and no stabilization is required in order to avoid spurious oscillations into the solution. 

Moreover, the spectral analysis is performed in the non-Hcrmitian case by re- 
ferring to the Hermitian and skew-Hermitian (HSS) decomposition of A n (a, /3) (that 
can be performed on any single elementary matrix related to Th by considering the 
standard assembling procedure). 

According to the definition, the HSS decomposition is given by 

Re(A n (aJ)) := Ma J) +^ (a J) = ^ + (2g) 
Hm(A n (aJ)) :=i A,(0 '^~^ ( °'^ =iIm(^)), (2-10) 
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where 

Re(*„(/3)) = \{* n 0) + *Z0)) = E n 0), (2.11) 

since by definition, the diffusion term 8„(a) is a Hermitian matrix and does not 
contribute to the skew-Hermitian part of A n (a,f3). Notice also that E n (f3) = if 
div(/3) = 0. More in general, the matrix Re(A n (a, (3)) is symmetric and positive def- 
inite whenever A m i n (0„(a)) > p(E n ((3)). Indeed, without the condition div(/3) > 0, 
the matrix E n (j3) does not have a definite sign in general: in fact, a specific analysis of 
the involved constants is required in order to guarantee the nonnegativity of the term 
E n ((3). Moreover, the Lemma below allows to obtain further information regarding 
such a spectral assumption, where || • H2, || • || 00 indicate both the usual vector norms 
and the induced matrix norms. 

Lemma 2.1. [10] Let {E n ((3)} be the matrix sequence defined according to (2.11). 
Under the assumptions in (2.2) 7 then we find \\E n (f3)\\2 < ||-En(/3)||oo < Ch 2 , with 
C absolute positive constant only depending on /3(x) and il. The claim holds both 
in the case in which the matrix elements in (2.7) are evaluated exactly and whenever 
a quadrature formula with error 0(h 2 ) is considered for approximating the involved 
integrals. 

Hereafter, we will denote by {A n (a, /?)}, n = n(h), the matrix sequence associated 
to a family of meshes {Th} 1 with decreasing finesse parameter h. As customary, the 
whole preconditioning analysis will refer to a matrix sequence instead to a single 
matrix, since the goal is to quantify the difficulty of the linear system resolution in 
relation to the accuracy of the chosen approximation scheme. 

3. Spectral analysis and clustering properties in the case of structured 
uniform meshes. The aim of this section is to analyze the spectral properties of 
the preconditioned matrix sequences {P~ 1 (a)(A n (a, /?))} in the case of some special 
domains Q, partitioned with structured uniform meshes, so that spectral tools derived 
from Toeplitz theory [6, 11] can be successfully applied. The applicative interest of 
the considered type of domains will be motivated in the short. 

Indeed, let / be a 2— variate Lebesgue integrable function defined over V = 
(— 7r,7r] 2 . By referring to the Fourier coefficients of this function /, called generating 
function, 

aj = i jj{s)e-' 1<j ' s> ds, i 2 = -1, 3 = (h,h) e Z 2 , 

with < j, s >= jisi + J2S2, one can build the sequence of Toeplitz matrices {T n (f)}. 
The matrix T n (f) e C nxn is said to be the Toeplitz matrix of order N = (JVi,7V 2 ), 
n = N1N2, generated by /. 

Now, the spectral properties of the matrix sequence {T n (f)} are completely un- 
derstood and characterized in terms of the underlying generating functions (see, e.g., 
[6] for more details). In the following, with respect to the latter claim, we will refer 
to the notion of equivalent generating functions. In such a respect, we claim that two 
nonnegative function / and g defined over a domain T> are equivalent if and only if 
/ = 0(g) and g = 0(f), where a — 0((3) means that there exists a pure positive 
constant c such that a < c(3 almost everywhere on V. 
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(a) (b) 

Fig. 3.f . Uniform structured mesh with equilateral or rectangular triangles. 

The main motivation of this paper lies in the analysis of the template case reported 
in Figure 3.1. a, where we consider a partition into equilateral triangles. It is self- 
evident that such a problem represents just an academic example. However, it is a 
fact that a professional mesh generator will locally produce a partitioning, which is 
"asymptotically" similar to the one reported in the figure. 

Our goal is to prove the PCG optimality with respect to the diffusion problem 
approximation, i.e., the number of iterations for reaching the solution, within a fixed 
accuracy, can be bounded from above by a constant independent of the dimension 
n = n(h). Some additional results about PGMRES convergence properties in the 
case of the convection-diffusion problem approximation are also reported, both in the 
case of the hexagonal domain with a structured uniform mesh as in Figure 3.1. a 
and in the case Q = (0, l) 2 with a structured uniform mesh as in Figure 3.1.b, so 
extending previous results proved in [10] with respect to the PHSS method [3]. 

The spectral analysis makes reference to the following definition. 

Definition 3.1. [23] Let {A n } be a sequence of Hermitian matrices of increasing 
dimensions n. The sequence {A n } is clustered at p in the eigenvalue sense if for any 
e > 0, #{« \Xi(A n ) tf: (p — e,p + e)} = o(n). The sequence {A n } is properly (or 
strongly) clustered at p if for any e > the number of the eigenvalues of A n not 
belonging to (p — e,p + e) can be bounded by a pure constant eventually depending on 
e, but not on n. In the case of non-Hermitian matrices, the same definition can be 
extended by considering the complex disk D(p,e), instead of the interval (p— e,p+s). 

3.1. Diffusion Equations. We start by considering the simpler diffusion prob- 
lem, i.e., the analysis concerns the Hermitian matrix sequences {A n (a}}, A n (a) = 

e n (a). 

Due to the very special choice of the domain O, the matrices A n (a) — 6„(a) 
arising in the case a = 1 with structured meshes as in Figure 3.1. a result to be a proper 
projection of Toeplitz matrices generated by the function /(si, S2) = \/3(6— 2 cos(si) — 
2 cos(s2) — 2 cos(si + S2))/3, (si,S2) E V = (— 7r,7r] 2 , related to a bigger parallelogram 
shaped domain fijv containing (see Figure 3.2), i.e., A n (l) = nTjv(/)n T , with N 
number of the internal nodes. Here, n € M™ xAr , n < N is a proper projection matrix 
simply cutting those rows (columns) referring to nodes belonging to £In, but not to 
f2. Thus, the matrix is full rank, i.e., rank(n) = n, and IH1 T = I n . 

In the same way, we can also consider the Toeplitz matrix with the same gener- 
ating function arising when considering a smaller parallelogram shaped domain Cl^ 
contained in tt, i.e., T^(f) = flA n (l)n T , with fl e M^ x ", N < n proper projection 
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Fig. 3.2. Uniform structured mesh on Q, Qm, 



matrix and such that rank(II) = N, and nn T = 1$. 

These embedding arguments allow to bound, according to the min-max principle 
(see [19] for more details on the use of this proof technique), the minimal eigenvalue 
of the matrices A n (l) as follows 

AminCM/)) < A min (A„(l)) < A min (T^(/)). (3.1) 

A first technical step in our spectral analysis concerns relationships between this 
generating function / and the more classical generating function /(si,S2) = 4 — 
2 cos(si) — 2 cos(s 2 ) arising in the case of FE approximations on a square f2 = (0, l) 2 
with Friedrichs-Keller meshes (see Figure 3.1.b), or standard FD discretizations. 

It can be easily observed that these two function are equivalent, in the sense of 
the previously reported definition, since we have 

^f<f<V3.f onV=(-7r,n} 2 . 

Thus, because T n {-) is a matrix- valued linear positive operator (LPO) for every n 
[13], the Toeplitz matrix sequences generated by this equivalent functions result to be 
spectrally equivalent, i.e., for any n 

^j-T n (f)<T n (f)<V3T n (f); (3-2) 

furthermore, due to the strict positivity of every T n (-) (see [13] for the precise defini- 
tion), since /— is not identically zero, we have that 

Tn(f)-^T n (f) 

is positive definite, and since \/3/ — / is not identically zero, we find that 

V3T n (f)-T n (f) 

is also positive definite. Here, we arc referring to the standard ordering relation 
between Hermitian matrices, i.e., the notation X > Y, with X and Y Hermitian 
matrices, means that X — Y is nonnegative definite. 

An interesting remark pertains to the fact that the function / is the most natural 
one from the FE point of view, since no contribution are lost owing to the gradient 
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orthogonality as instead in the case related to /; nevertheless its relationships with the 
function / can be fully exploited in performing the spectral analysis. More precisely, 
from (3.1) and (3.2) and taking into account (see e.g. [11]) that 

A mi „(T„(/)) = 8 sin 2 f^h) , h = — |— , 
V 2 / n + 1 

we deduce 

A min (A„(l)) > A min (TAr(/)) > \ X min (T N (f)) = ^-8 sin 2 (|^) ~ h 2 N . (3.3) 
Following the very same reasoning, we also find that 

A min (A„(l)) < A min (7>(/)) < \^A min (T JV (/)) = 873sin 2 (|/^) ~ 4. (3.4) 

Since, by the embedding argument, both N and N are asymptotic to n, it follows that 
A m in(Ai(l)) ~ ^ 2 - Finally, following the same analysis for the maximal eigenvalue, 
we find 

/3 

^- A max (T JV (/)) < A max (A„(l)) < V3 A max (T^(/)), 

where, by [6], we know that 

lim A max (T n (/)) = max / = 8, 

n— >oo T> 

and hence the spectral condition number of A n (l) grows as h~ 2 i.e. 

K 2 (A n (l))~h- 2 , A n (l) = 6„(l), (3.5) 

where the constant hidden in the previous relation is mild and can be easily estimated. 

It is worth stressing that the same matrix II can also be considered in a more 
general setting. In fact, the matrix sequence {A n (a)} can also be defined as {A n (a)} = 
{IIAjv(a) II T }, since again each internal node in £1 is a vertex of the same constant 
number of triangles. Therefore, by referring to projection arguments, the spectral 
analysis can be equivalently performed both on the matrix sequence {A n (a)} and 
{A N (a)}. 

No matter about this choice, we make use of a second technical step, which is 
based on standard Taylor's expansions. 

Lemma 3.2. Let a e C 2 (fl), with a(x) > ao > and hexagonal domain 
partitioned as in Figure 3.1. a. For any p such that (©n(a)) s ,s-p 7^ there exists a 
proper (x*,y*) such that the Taylor's expansions centered at a proper (x*,y*) have 
the form 

(®n(a))s,s- P = a(x* p , y*)(A n (l)) s _ p + h 2 D p + o(h 2 ), (3.6) 

(e„(o)),,. = a(x;, y ;)(A n (i)) s , s + hB p + h 2 c p + o(h 2 ), (3.7) 

(6„(a)) s _ p , s _ p - a(x* p , y ;)(A n (l)) s ,s - hB p + h 2 C p + o(h 2 ), (3.8) 

where B P ,C P and D p are constants independent of h. 

Proof. In the cases at hand, the validity of this claim is just a direct check: the key 
point relies in the symmetry properties induced by the structured uniform nature of 
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Fig. 3.3. Centers of Taylor's expansion. 

the considered meshes, both in the case of f2 and fijy Hereafter, we give some detail 
with respect to the case of the hexagonal domain partitioned as in Figure 3.1. a, see 
[17] for the proof in the case Q. = (0, l) 2 . Thus, let's consider the Taylor's expansion 
centered at a proper (x*,y*). By calling a* = a(x*,y*) and similarly denoting the 
derivatives of the diffusion coefficient, it holds that 

/ a(x,y)= a* + (x - x*)a* x + (y - y*)a* y 

JK JK 

+ \{x - x *?< x + - (y - y*) 2 a* yy + (x - x*)(y - y*)a* xy + o(t 2 ) 
= \ K \(a* +{x- x*)\ bK a* x + (y- y*)\ bK a* y ) 

/I i 
r 2^ - x *f<x + ^(y - y*) 2a * yy + ( x ~ x *)(y - y*)a* xy + o(r 2 ) 

= \ K \( a * +(x- x*)\ bK a* x + {y- y*)\b K a* y 

+\( x ~ x % K <x + \iy- y% K a ly + ( x - x *)(y - y*)\b K < y ) + o(t 2 ) 

with t = \\[x — x'y — y*] T \\ and denoting the barycenter of the triangle K. By 

referring to the assembling procedure, we will choose as center of the related Taylor's 

expansion in the case of Figure 3.3.b-3.3.d the point marked by *. The symmetric 

position of the involved triangle barycenters (marked by o), with respect to * allows 

to claim (3.6). The same holds true for (3.7) by considering each pertinent point * in 

Figure 3. 3. a and, analogously, for (3.8). □ 

Lemma 3.3. [17] Let's considering the following notation of scaled matrix X* = 
_i _ i 

D n 2 (a)XD n 2 (a) of a given matrix X. Under the assumptions of Lemma 3.2 the 
following representation holds true 

e* n (a) = ^- 1 / 2 ( a )6„(a) J D- 1 / 2 ( a ) = 0„(1) + h 2 F n (a) + o{h 2 )G n {a) (3.9) 

where F n and G n are uniformly bounded in spectral norm and have the same pattern 
as 0*(a). 
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We are now ready to prove the PCG optimality. 

Theorem 3.4. Let {A n (a)} and {P n (a)} be the Hermitian positive definite ma- 
trix sequences defined according to (2.5) and (2.8) in the case of the hexagonal do- 
main partitioned as in Figure 3.1. a. Under the assumptions in (2.2), the sequence 
{P~ (a)A n (a)} is properly clustered at 1. Moreover, for any n all the eigenvalues of 
P^ 1 (a)A n (a) belong to an interval [d, D] well separated from zero [Spectral equiva- 
lence property]. 

Proof. Since A n (a) = nA N (a)Il T , with II e R nxN such that rank(II) = n, and 
II II T = I n , we have 

p-\a)A n {a) = [L»|(a)A„(l,0)£>|(a)]- 1 A„(a) 

= [n £>* ( a ) n T iL4jv(i) n T n £>|(a) n T ]- x n A N (a) n T 

= [nn T n Dl(a)A N (i)D\{a) n T n n T ]- x n A N (a) n T 

= [n£)J(a)^(i)£>| r (a)n T ]- 1 n^ jV (a)n T 
= pp^ 1 (o)n r ]- 1 n^(a)n r . 

Since II is full rank, it is evident that the spectral behavior of 

[np A r( a )n T ]- 1 nA w (a)n T 

is in principle better than the one of P^ 1 (a)Aisf(a), to which we can address the 
spectral analysis. Thus, the proof technique refers to [17, 19] and the very key step 
is given by the relations outlined in (3.3) and (3.4). □ 

It is worth stressing that the previous claims can also be proved by considering 
the sequence {P~ x (a)A n (a)} , instead of {P^ 1 (a)Ajv(a)}, simply by directly referring 
to the quoted asymptotical expansions. The interest may concern the analysis of the 
same uniform mesh on more general domains. 

Moreover, the same spectral properties has been proved in the case of uniform 
structured meshes as in Figure 3.1.b in [17]. 

3.2. Convection-Diffusion Equations. The natural extension of the claim in 
Theorem 3.4, in the case of the matrix sequence {Re(A n (a, /?))} with Re(A n (a, /?)) ^ 
0„(a), can be proved under the additional assumptions of Lemma 2.1, in perfect 
agreement with Theorem 5.3 in [10]. 

Theorem 3.5. Let {Re(A n (a, (3))} and {P n (a)} be the Hermitian positive defi- 
nite matrix sequences defined according to (2.9) and (2.8) in the case of the hexagonal 
domain partitioned as in Figure 3.1. a. Under the assumptions in (2.2), the sequence 
{P~ 1 (a)Re(A n (a, /?))} is properly clustered at 1. Moreover, for any n all the eigen- 
values of P~ 1 (a)Re(A n (a, (3)) belong to an interval [d,D] well separated from zero 
[Spectral equivalence property] . 

The claim holds both in the case in which the matrix elements in (2.6) and (2.7) 
are evaluated exactly and whenever a quadrature formula with error 0(h 2 ) is consid- 
ered to approximate the involved integrals. 

Proof. The proof can be done verbatim as in [10], since the key point is proved 
by referring to relations in (3.3) and (3.4). □ 

Theorem 3.6. Let {lm(A n (a, (3))} and {P n (a)} be the Hermitian matrix se- 
quences defined according to (2.10) and (2.8) in the case of the hexagonal domain 
partitioned as in Figure 3.1. a. 

Under the assumptions in (2.2), the sequence {P~ 1 (a)lm(A n (a, /?))} is spectrally 
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bounded and properly clustered at with respect to the eigenvalues. The claim holds 
both in the case in which the matrix elements in (2.7) are evaluated exactly and when- 
ever a quadrature formula with error 0{h 2 ) is considered to approximate the involved 
integrals. 

Proof. This result has been proved in [10] with respect to Friedrichs-Keller trian- 
gulations, but it can easily be extended to the case of matrix sequences arising in the 
case of a FE partitioning as in Figure 3.1. a by using the same arguments. □ 

On the basis of these two splittcd spectral results, we can easily obtain the spectral 
description of the whole preconditioned matrix sequence {P~ 1 (a)A n (a, /?)}, according 
to the theorem below. The proof technique refers to an analogous result proved in [4] 
with respect to FD discretizations of convection-diffusion equations. 

Theorem 3.7. Let {(A n (a, /?))} and {P n (a)} be the matrix sequences defined 
according to (2.5) and (2.8), both in the case of the hexagonal domain Q with a struc- 
tured uniform mesh as in Figure 3.1. a and in the case = (0, l) 2 with a structured 
uniform mesh as in Figure S.l.b 

Under the assumptions in (2.2), the sequence {P~ 1 (a)A n (a, ft)} is properly clus- 
tered at 1 € C + with respect to the eigenvalues. In addition, these eigenvalues all 
belong to a uniformly bounded rectangle with positive real part, well separated from 
zero. 

The claim holds both in the case in which the matrix elements in (2.6) and (2.7) 
are evaluated exactly and whenever a quadrature formula with error 0(h 2 ) is consid- 
ered to approximate the involved integrals. 

Proof. A localization result for the eigenvalues of the sequence {P~ 1 (a) (A n (a,$))} 
can be stated simply by referring to the properties of the field of values of the pre- 
conditioned matrix: in fact, for any i 

XiiP-^^A^aJ) ET = F (p- 1 -{a)A n {aJ)p- 1 -{a)) 



where 



jr. J, e C : , _ .%(4M)). +i /l.W. : ai» i x £ c „ m 



x H P n (a)x x H P n (a)x 



Since for any i 



X i ( P -\a)Re(A n (aJ)) e jz € C : z = ^M^ , x e C"\{0}| , 

K{ P -\a)l m {A n {aJ)) e jz € C : z = ^^f^ , * e C«\{0}| . 

the claimed localization result is obtained as a consequence of those previously proved 
in Theorems 3.5 and 3.6. Moreover, the same localization results together with the 
claimed clustering properties allow to apply Theorem 4.3 in [4] and the proof is con- 
cluded. □ 

The subsequent result, giving estimates on the condition number of the eigenvec- 
tor matrix of the preconditioned structure, is of interest in the study of the convergence 
speed of the GMRES. In particular a logarithmic growth in number of iteration is as- 
sociated to a polynomial bound in the spectral condition number: such a bound is 
precisely given below. 
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Theorem 3.8. Let a(x) and /?(x) constant functions and let {(A n (a, /?))} and 
{P„(a)} 6e the matrix sequences defined as in Theorem 3.7 
Then, P„(a) _1 A n (a, /3) can 6e diagonalized as 

P n (a)- 1 A n (aJ)=V n D n V- 1 

where the matrix of the eigenvectors V n can be chosen such that K2(V n ) ~ 

Proof. Under the quoted assumption it holds that P n {a) = a© n (l) and is 
skew-symmetric. Thus, we have 

P n {a)- l A n (aJ) = / n + a- 1 6- 1 (l)*„(/3) 

= / n + e; l (i)w n e|(i), 

_1 _i -> _1 _1 

with W n = a~i 0„ 2 (l)W n (/3)0„ 2 (l)a~ 5 skew-symmetric matrix. Therefore, the 

matrix J„ + W n is a normal normal, so that /„ + W n — Q n DnQn , where D n is 

diagonal and Q n is unitary. Consequently, 

P n (a)- x A n {aJ) = 0^'(l)(J n + VK„)e|(l) 
= e~*(l)Q n £> n Q*e*(l) 
= KA^" 1 

_ i_ 

with V n = O n 2 (l)Q n eigenvector matrix. Due to the fact that Q n is unitary, we have 

K 2 (V n ) = K 2 {Qn^{l)Q n ) = K 2 {@~^{1)) and finally K 2 (V n ) - h' 1 by invoking the 
key relation (3.5). □ 

4. Numerical tests. The section is divided into two parts. In the first we briefly 
discuss the complexity features of our preconditioning proposals. In the second part we 
report and critically discuss few numerical experiments in which the meshes are both 
structured and unstructured, while the regularity features required by the theoretical 
analysis are somehow relaxed. 

4.1. Complexity issues. First of all, we report some remarks about the com- 
putational costs of the proposed iterative procedure. 

The main idea is that such a technique is easily applicable. In fact, a Krylov 
method is considered, so simply requiring a matrix vector routine for sparse matrices 
and a solver for the chosen preconditioner. 

Therefore, since the preconditioner is defined as P n (a) = D 1 J 2 (a)A n (l,0)Dn 2 (a), 
where D n {a) = diag(A„(a, 0))diag _1 ( J 4„(l, 0)), the solution of the linear system in 
(2.4) with matrix A n (a,f3) is reduced to computations involving diagonals and the 
matrix A n (l,0) (A n (a) and A n (l,0) for the diffusion problem, respectively). 

As well known, whenever the domain is partitioned by considering a uniform 
structured mesh this latter task can be efficiently performed by means of fast Poisson 
solvers, among which we can list those based on the cyclic reduction idea (see e.g. 
[7, 8, 21]) and several specialized algebraic or geometric multigrid methods (see e.g. [9, 
22, 14]). In addition, the latter can be efficiently considered also in more general mesh 
settings. The underlying idea is that the main effort in devising efficient algorithms 
must be devoted only to this simpler problem with constant coefficient, instead of the 
general one. 
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Now, the effectiveness of the proposed method is measured by referring to the 
optimality definition below. 

Definition 4.1. [1] Let {A m x m = b m } be a given sequence of linear systems of 
increasing dimensions. An iterative method is optimal if 

1. the arithmetic cost of each iteration is at most proportional to the complexity 
of a matrix-vector product with matrix A m , 

2. the number of iterations for reaching the solution within a fixed accuracy can 
be bounded from above by a constant independent ofm. 

In such a respect Theorem 3.4 proves the optimality of the PCG method: the 
iterations number for reaching the solution within a fixed accuracy can be bounded 
from above by a constant independent of the dimension n = n(h) and the arithmetic 
cost of each iteration is at most proportional to the complexity of a matrix- vector 
product with matrix A n (a, 0). 

Moreover, Theorem 3.7 proves that all the eigenvalues of the preconditioned ma- 
trix belong in a complex rectangle {z e C : Rc(z) e [d, D], Im(z) e [— d, d]}, with 
d, D > 0, d > independent of the dimension n = n(h). It is worth stressing that 
the existence of a proper eigenvalue cluster and the aforementioned localization re- 
sults in the preconditioned spectrum can be very important for fast convergence of 
preconditioned GMRES iterations (see, e.g., [5]). 

Finally, we want to give notice that the PCG/PGMRES numerical performances 
do not get worse in the case of unstructured meshes, despite the lack of a rigorous 
proof. 

4.2. Numerical results. Before analyzing in detail selected numerical results, 
we wish to give technical information on the performed numerical experiments. We 
apply the PCG or the PGMRES method, in the symmetric and non-symmetric case, 
respectively, with the preconditioning strategy described in section 2, to FE approx- 
imations of the problem (1.1). Whenever required, the involved integrals have been 
approximated by means of the barycentric quadrature rule (the approximation by 
means of the nodal quadrature rule gives rise to similar results, indeed both are exact 
when applied to linear functions). 

The domains of integration fi are those reported in Figure 3.1 and we assume 
Dirichlet boundary conditions. All the reported numerical experiments are performed 
in Matlab, by employing the available peg and gmres library functions; the iterative 
solvers start with zero initial guess and the stopping criterion 1 1 r-fe 1 1 2 < 1 — T 1 1 r"o 1 1 2 is 
considered. The case of unstructured meshes is also discussed and compared, together 
with various types of regularity in the diffusion coefficient. 

In fact, we consider the case of a coefficient function satisfying the assumptions 
(2.2). More precisely, the second columns in Table 4.1 report the number of iterations 
required to achieve the convergence for increasing values of the coefficient matrix size 
n = n{h) when considering the FE approximation with structured uniform meshes as 
in Figure 3.1 and with template function a(x,y) — a\(x,y) — exp(x + y), (3(x,y) = 
[x y] T .The subsequent meshes are obtained by means of a progressive refinement 
procedure, consisting in adding new nodes corresponding to the middle point of each 
edge. The numerical experiments plainly confirm the previous theoretical analysis 
in section 3: the convergence behavior does not depend on the coefficient matrix 
dimension n = n(h). 

As anticipated, despite the lack of corresponding theoretical results, we want to 
test the convergence behavior also in the case in which the regularity assumption 
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Table 4.1 

Number of PCG and PGMRES iterations - structured meshes as in Fig. 3.1. a, Q hexagonal 
domain. 



PCG 


PGMRES, ^(x,j/) = [x 




n 


a x {x,y) 


a 2 (x,y) 


a-i(x,y) 


n 


ai(x,y) 


a 2 (x,y) 


03(2,2/) 


37 


3 


4 


5 


37 


4 


5 


5 


169 


3 


4 


5 


169 


4 


5 


5 


721 


3 


4 


4 


721 


4 


5 


5 


2977 


3 


4 


4 


2977 


4 


5 


5 


12097 


3 


4 


4 


12097 


4 


5 


5 


48769 


3 


4 


4 


48769 


4 


5 


5 



on a(x, y) in Theorems 3.4 and 3.7 are not satisfied. The analysis is motivated by 
favorable known numerical results in the case of FD approximations (see, for instance, 
[18, 19, 3]) or FE approximation with only the diffusion term [17]. More precisely, we 
consider as template the C 1 function a{x, y) = a^ix, y) = e x+ \ y ~ y °\ 3/2 , the C° function 
a(x,y) = a 3 {x,y) = e^+lv-Kol, wit h y = vT3)/4 or y Q = 1/2. 

The number of required iterations is listed in the remaining columns in Table 4.1. 
Notice that also in these cases with a C 1 or C° diffusion function, the iteration count 
does not depend on the coefficient matrix dimension n = n(h). The same results 
are reported in Table 4.2 with respect to structured uniform meshes on the domain 
f2 = (0,l) 2 . 

Furthermore, we want to test our proposal in the case of some unstructured 
meshes generated by triangle [20] with a progressive refinement procedure. The first 
meshes in the considered mesh sequences are reported in Figures 4.1 and 4.2, respec- 
tively. 

Tables 4.3 and 4.4 report the number of required iterations in the case of the 
previous template functions. Negligible differences in the iteration trends are observed 
for increasing dimensions n. All these remarks arc in perfect agreement with the 
outliers analysis of the matrices P~ 1 {a)A n {a 1 /3), with respect to a cluster at 1 e C + 
(see some examples in Figure 4.3). 

To conclude the section we take into consideration the more realistic setting in 
which the meshes are generated by a specialized automatic procedure. According to 
this point of view, we have applied a spectral approximation of the matrix sequence 
{A n (a)} in terms of product of low-cost matrix structures, i.e., A n (l) and D n (a), that 
carry the "structural" content of the variational problem and the specific "informa- 
tive" content contained in the diffusion coefficient function a, respectively. However, 
the matrix A n (l) can be still not easy to handle. Therefore, we go beyond in such 
idea by defining a preconditioner P n (a) in which the matrix A n (l) related to the 
unstructured mesh is replaced by a suitable projection of the Toeplitz matrix Tjv(/), 
f(si, 82) = V3(6 - 2 cos(si) - 2 cos(s 2 ) - 2 cos(si + S2))/3, (si, S2) eP = (— 7r, it] 2 
(see Section 3.1). 

The numerical results are reported in Table 4.5. As expected, the number of 
iterations is a constant with respect to the dimension when the preconditioner P„(a) 
is applied. Nevertheless, for n large enough, the same seems to be true also for the 
preconditioner P n {a). Indeed, for increasing dimensions n, the unstructured parti- 
tioning is more and more similar to the one given by equilateral triangles sketched in 
Fig. 3.1.a. 

A theoretical ground supporting these observation is still missing and would be 
worth in our opinion to be studied and developed. 
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Table 4.2 

Number of PCG and PGMRES iterations - structured meshes as in Fig. 3.1.b, fl = (0, l) 2 . 



PCG 


PGMRES, ${x,y) = [x y] T 


n 


ai(x,y) 


a 2 (x,y) 


a.3{x, y) 


n 


ai(x,y) 


a 2 (x,y) 


0.3 (x, y) 


81 


3 


4 


4 


81 


4 


5 


5 


361 


3 


4 


5 


3611 


4 


5 


5 


1521 


3 


4 


5 


1521 


4 


5 


5 


6241 


3 


4 


5 


6241 


4 


5 


5 


25281 


3 


4 


5 


25281 


4 


5 


5 


128881 


3 


4 


5 


128881 


4 


5 


5 



Table 4.3 

Number of PCG and PGMRES iterations - unstructured meshes on the hexagonal domain Q. 



PCG 


PGMRES, /3(x,j0 = [x 




n 


a 1 (x,y) 


a 2 (x,y) 


a,3{x,y) 


n 


ai(x,y) 


a 2 (x,y) 




28 


5 


5 


5 


28 


4 


5 


5 


73 


4 


4 


5 


73 


4 


5 


5 


265 


4 


4 


5 


265 


4 


5 


5 


1175 


4 


4 


5 


1175 


4 


5 


5 


4732 


4 


4 


5 


4732 


4 


5 


5 


19288 


4 


4 


5 


19288 


4 


5 


5 


76110 


4 


4 


4 


76110 


4 


5 


5 



Table 4.4 

Number of PCG and PGMRES iterations - unstructured meshes on the domain Q = (0, l) 2 . 



PCG 


PGMRES, f3(x,y) = [x 




n 


ai(x,y) 




a,3{x,y) 


n 


ai(x,y) 


a 2 (x,y) 


a 3 (x, J/) 


24 


4 


5 


5 


24 


4 


5 


5 


109 


4 


5 


5 


109 


4 


5 


5 


465 


4 


5 


5 


465 


4 


5 


5 


1921 


4 


5 


5 


1921 


4 


5 


5 


7809 


4 


5 


5 


7809 


4 


5 


5 


31489 


4 


5 


5 


31489 


4 


5 


5 



Table 4.5 

Number of PCG and PGMRES iterations - structured and unstructured meshes on the hexag- 
onal domain Q, a\{x,y), ft(x,y) = [x y] T . 



PCG 




PGMRES 


n 


P 


P 


n 


P 


P 


37 


4 


8 




37 


4 


8 


169 


4 


10 




169 


4 


9 


721 


4 


11 




721 


4 


9 


2977 


4 


11 




2977 


4 


10 


12095 


4 


13 




12095 


4 


11 


48769 


4 


16 




48769 


4 


13 
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Fig. 4.1. Unstructured Meshes on the hexagonal domain U. 




1G 



A. RUSSO, S. SERRA CAPIZZANO, AND C. TABLINO POSSIO 




Fig. 4.3. Outlier analysis - structured, perturbed and unstructured meshes on the hexagonal 
domain f2, ai(x,y), a,2(x,y), as{x,y), and j3{x,y) = [x y] T . 



5. Perspectives and future works. As emphasized in the introduction it is 
clear that the problem in Figure 3.1. a is just an academic example, due to the perfect 
structure made by equi-lateral triangles. However, it is a fact that a professional mesh 
generator will produce a partitioning, which "asymptotically", that is a for a mesh 
fine enough, tends to the one in Figure 3.1. a. 

The latter fact has a practical important counterpart, since the academic precon- 
ditioner P n (a) is optimal for the real case with nonconstant coefficients and with the 
partitioning in B). A theoretical ground supporting these observations is still missing 
and would be worth in our opinion to be studied and developed in three directions: a) 
giving a formal notion of convergence of a partitioning to a structured one, b) proving 
spectral and convergence results in the case of an asymptotically structured partition- 
ing, c) extending the spectral analysis in the case of weak regularity assumptions. 

Other possible developments include the case of higher order finite element spaces: 
it would be intriguing to find an expression of the underlying Toeplitz symbol as a 
function of the different parameters in the considered finite element space (degrees of 
freedom, polynomial degree, geometry of the mesh), and this could be done uniformly 
in d dimension, i.e. for equation (1.1) with C R d , d>2. 
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